% 说明:
% 1 小波变换使用 a'trous 快速小波提升算法(即教材P96 P188提及的 多孔算法) 见参考文献 《二维A Tuous算法图像边缘检测研究》;
% 2 二维二进小波变换快速算法的推导过程 见教材P195  式(7-42)是推导结论;

% 目前已解决问题:
% 1 二维二进小波滤波器的计算构造(实质是求两尺度方程系数 目前没有成熟的计算方法 
%   仅能查到结果并且有纰漏 选取教材提供的两组数据作为解决方案)
%     1.1  非正交二次样条小波
%     1.2  正交二次样条小波(近似解与修正解)
% 2 在确定阈值T的情况下边缘检测
%     2.1  根据滤波器系数进行小波变换求像素点的梯度值和梯度角
%     2.2  非极大值抑制
%     2.3  阈值化
%     2.4  边缘连接

% 目前遗留问题:
% 3 阈值的自适应算法 或者通过GUI作为手工输入参数处理
% 4 多个尺度边缘下结果的筛选处理
% 5 工程问题: 图形界面中应该设置哪些参数及实现;

clear all; 
X=imread('Lena.jpg');
I=rgb2gray(X);
figure;
imshow(I); 
[N,M] = size(I); 

%以下h g delta三个变量的意义: 分别对应于教材P197(7-42) 的 h g l滤波器系数

%1 非正交二次样条小波系数的构造 参考 教材P182例7.1和P193例7.4 需要阈值化得到边缘
% h = [0.125,0.375,0.375,0.125]; 
% g = [-1.0,1.0]; 
% delta = [0.015625,0.09375,0.234375,1.3125,0.234375,0.09375,0.015625];

%2 正交二次样条小波系数的构造 参考 教材P183例7.2和P194例7.5 需要阈值化得到边缘
h = [0.125,0.375,0.375,0.125]; 
%g = [0.0061,0.0869,0.5798,-0.5798,-0.0869,-0.0061];%近似解 
g = [-0.9324,0.6736,0.2253,0.0335];%修正解 from PPT
delta = [0.00003,0.00727,0.03118,0.06623,0.79113,0.06623,0.03118,0.00727,0.00003];

%参数K意义:多尺度边缘检测中的尺度个数
K = 4; 
%tan_angle = zeros(N, M); %角度矩阵(此处采用角度，也可直接用正切来比较)

%四维矩阵  实际上只使用到其中的后两维存储二维小波变换的结果
a(1:N,1:M,1,1:K+1) = 0; 
dx(1:N,1:M,1,1:K+1) = 0; 
dy(1:N,1:M,1,1:K+1) = 0; 
d(1:N,1:M,1,1:K+1) = 0; 
tan_angle(1:N,1:M,1,1:K+1) = 0;


%边缘化第一步 对I进行处理
a(:,:,1,1) = conv2(h,h,I,'same'); 
dx(:,:,1,1) = conv2(delta,g,I,'same'); 
dy(:,:,1,1) = conv2(g,delta,I,'same'); 

x = dx(:,:,1,1); 
y = dy(:,:,1,1); 

d(:,:,1,1) = sqrt(x.^2+y.^2); %计算1尺度下的梯度模值
dt=d(:,:,1,1);%dt at为中间变量
tan_angle(:,:,1,1) = atan(x/y);%计算1尺度下的梯度角
at=tan_angle(:,:,1,1);



%沿着梯度方向进行 非极大值抑制
PT=dt;%dt在非极大值抑制的过程中会发生改变故将原始dt存于PT 然后基于PT对dt进行非极大值抑制
for I = 2:N-1
    for J = 2:M-1
        if at(I, J) >= -pi/8 && at(I, J) < pi/8 && ( PT(I, J)...
                <=  PT(I+1, J) || PT(I, J) <=  PT(I-1, J))
            dt(I, J) = 0;
        end
        if at(I, J) >= pi/8 && at(I, J) < 3*pi/8 && ( PT(I, J)...
                <=  PT(I-1, J+1) ||  PT(I, J) <=  PT(I+1, J-1))
             dt(I, J) = 0;
        end
        if at(I, J) >= -3*pi/8 && at(I, J) < -pi/8 && ( PT(I, J)...
                <=  PT(I-1, J-1) ||  PT(I, J) <=  PT(I+1, J+1))
            dt(I, J) = 0;
        end
        if (at(I, J) < -3*pi/8 || at(I, J) >=  3*pi/8) && ( PT(I, J)...
                <=  PT(I-1, J) ||  PT(I, J) <=  PT(I+1, J))
            dt(I, J) = 0;
        end
    end
end

%阈值化
P=zeros(N,M);
for i=1:N
    for j=1:M
        if dt(i,j)>25
            P(i,j)=255;
        end
    end
end

%边缘连接
PT=P;
for I = 2:N-1
    for J = 2:M-1
        if at(I, J) >= -pi/8 && at(I, J) < pi/8 && ( PT(I, J)...
                == 255)
            P(I+1, J)=255;
            P(I-1, J)=255;
            dt(I, J) = 0;
        end
        if at(I, J) >= pi/8 && at(I, J) < 3*pi/8 && ( PT(I, J)...
                == 255)
            P(I-1, J+1)=255;
            P(I+1, J-1)=255;
        end
        if at(I, J) >= -3*pi/8 && at(I, J) < -pi/8 && ( PT(I, J)...
                ==  255)
            P(I-1, J-1)=255;
            P(I+1, J+1)=255;
        end
        if (at(I, J) < -3*pi/8 || at(I, J) >=  3*pi/8) && ( PT(I, J)...
                ==  255)
            P(I-1, J)=255;
            P(I+1, J)=255;
        end
    end
end

figure;
P=uint8(P);
imshow(P);

lh = length(h); 
lg = length(g); 

%接着进行J个尺度下的边缘检测 每次对上一次的低频a(:,:,1,j)进行处理
for j = 1:K
    lhj = 2^j*(lh-1)+1; 
    lgj = 2^j*(lg-1)+1; 
    hj(1:lhj)=0; 
    gj(1:lgj)=0; 
    for n = 1:lh 
        hj(2^j*(n-1)+1)=h(n); 
    end

    for n = 1:lg
        gj(2^j*(n-1)+1)=g(n); 
    end

    a(:,:,1,j+1) = conv2(hj,hj,a(:,:,1,j),'same');
    dx(:,:,1,j+1) = conv2(delta,gj,a(:,:,1,j),'same'); 
    dy(:,:,1,j+1) = conv2(gj,delta,a(:,:,1,j),'same'); 

    x = dx(:,:,1,j+1); 
    y = dy(:,:,1,j+1); 
    
    dj(:,:,1,j+1) = sqrt(x.^2+y.^2); %计算j+1尺度下的梯度模值
    dt=dj(:,:,1,j+1);
    tan_angle(:,:,1,j+1) = atan(x/y);%计算j+1尺度下的梯度角
    at=tan_angle(:,:,1,j+1);
    
    %沿着梯度方向进行非极大值抑制
    PT=dt;
    for I = 2:N-1
        for J = 2:M-1
            if at(I, J) >= -pi/8 && at(I, J) < pi/8 && ( PT(I, J)...
                    <=  PT(I+1, J) || PT(I, J) <=  PT(I-1, J))
                dt(I, J) = 0;
            end
            if at(I, J) >= pi/8 && at(I, J) < 3*pi/8 && ( PT(I, J)...
                    <=  PT(I-1, J+1) ||  PT(I, J) <=  PT(I+1, J-1))
                dt(I, J) = 0;
            end
            if at(I, J) >= -3*pi/8 && at(I, J) < -pi/8 && ( PT(I, J)...
                    <=  PT(I-1, J-1) ||  PT(I, J) <=  PT(I+1, J+1))
                dt(I, J) = 0;
            end
            if (at(I, J) < -3*pi/8 || at(I, J) >=  3*pi/8) && ( PT(I, J)...
                    <=  PT(I-1, J) ||  PT(I, J) <=  PT(I+1, J))
                dt(I, J) = 0;
            end
        end
    end
    %阈值化
    P=zeros(N,M);
    for i=1:N
        for j=1:M
            if dt(i,j)>25
                P(i,j)=255;
            end
        end
    end
    %边缘连接
    PT=P;
    for I = 2:N-1
        for J = 2:M-1
            if at(I, J) >= -pi/8 && at(I, J) < pi/8 && ( PT(I, J)...
                    == 255)
                P(I+1, J)=255;
                P(I-1, J)=255;
                dt(I, J) = 0;
            end
            if at(I, J) >= pi/8 && at(I, J) < 3*pi/8 && ( PT(I, J)...
                    == 255)
                P(I-1, J+1)=255;
                P(I+1, J-1)=255;
            end
            if at(I, J) >= -3*pi/8 && at(I, J) < -pi/8 && ( PT(I, J)...
                    ==  255)
                P(I-1, J-1)=255;
                P(I+1, J+1)=255;
            end
            if (at(I, J) < -3*pi/8 || at(I, J) >=  3*pi/8) && ( PT(I, J)...
                    ==  255)
                P(I-1, J)=255;
                P(I+1, J)=255;
            end
        end
    end
    
    
    figure;
    P=uint8(P);
    imshow(P);
end
